    def get_asa_old(self, query_ents =None, expand_symmetry =False, expand_crystal =False, probe =1.4, points =100, bucket_size =10):
        """ Calculates the atomic solvent-accessible surface area (ASA).
        """
        const = 4.0 * pi / points       # to normalize area per point
        points = sphere_points(points)  # point coordinates randomly distributed on the unit sphere.

        # The default is to
        lattice_res = einput(self, 'R')
        lattice_res._strip_children('H_HOH', 'eq', 'name')

        qents = einput((query_ents or lattice_res), 'A')
        qents.dispatch('set_radius')
        qents_radii = array(qents._data_children('radius')) + probe
        qents_idx_to_id = dict(enumerate(qents._data_children('get_full_id', forgiving =False, function =True)))
        qcoords = array(qents._data_children('coords', forgiving =False)) # keep them orthogonal

        lents = einput(lattice_res, 'A')
        lents.dispatch('set_radius')
        lents_radii = array(lents._data_children('radius')) + probe
        lents_idx_to_id = dict(enumerate(lents._data_children('get_full_id', forgiving =False, function =True)))
        lcoords = array(lents._data_children('coords', forgiving =False))

        if expand_symmetry:
            lcoords = self.coords_to_symmetry(lcoords, self.fmx, self.omx, self.mxs, self.raw)
        else:
            lcoords = array([lcoords])
        if expand_crystal:
            lcoords = self.coords_to_crystal(lcoords, self.fmx, self.omx, self.n)
        else:
            lcoords = array([lcoords])

        shape = lcoords.shape
        lcoords = lcoords.reshape((shape[0] * shape[1] * shape[2], shape[3]))
        kdt = kdtree(lcoords, bucket_size)

        result = {}
        for idx, qcoord in enumerate(qcoords):
            qradius = qents_radii[idx]
            #print qradius
            lradius_max = 2.0 + probe
            #print lradius_max
            search_limit = (qradius + lradius_max) * (qradius + lradius_max)
            #print qcoord
            #print search_limit
            ridxs, rdists =  kdt.rn(qcoord, search_limit)
            #print len(ridxs)
            ridxs = ridxs[rdists > 0]

            n_acc_points = 0
            for point in points:
                is_accessible = True
                real_point = point * qradius + qcoord           # expand and shift position
                #print '\t', point, real_point
                for ridx in ridxs:
                    entity_idx = ridx % shape[2]
                    lradius = lents_radii[entity_idx]
                    #print '        ',entity_idx, lradius,lcoords[ridx],
                    diff = real_point - lcoords[ridx]
                    #print diff,
                    distance_sq = diff * diff
                    #print distance_sq
                    if distance_sq.sum() < lradius * lradius:
                        is_accessible = False
                        break
                if is_accessible:
                    n_acc_points +=1
            result[qents_idx_to_id[idx]] = const * n_acc_points * qradius * qradius
        return result

    def get_contacts_old(self, query_ents =None, lattice_ents =None, level ='A', expand_symmetry =True,\
                        expand_crystal =True, search_limit =5.0, mode ='diff_sym',\
                        bucket_size =10):
        """ Search contacts between an query and target entity sets.
            mode:   diff_asu - require diff. asu
                    diff_sym - require diff. symmetry operator if the same asu (default)
                    diff_chain - require diff. chain if the same asu and symmetry operator

        """
        query_ents = (query_ents or self[(0,)])
        lattice_ents = (lattice_ents or self[(0,)])
        search_limit = search_limit * search_limit

        lents = einput(lattice_ents, level)
        lents_idx_to_id = dict(enumerate(lents._data_children('get_full_id', forgiving =False, function =True)))
        lcoords = array(lents._data_children('coords', forgiving =False))

        qents = einput(query_ents, level)
        qents_idx_to_id = dict(enumerate(qents._data_children('get_full_id', forgiving =False, function =True)))
        qcoords = array(qents._data_children('coords', forgiving =False)) # keep them orthogonal

        if expand_symmetry:
            lcoords = self.coords_to_symmetry(lcoords, self.fmx, self.omx, self.mxs, self.raw)
        else:
            lcoords = array([lcoords])
        if expand_crystal:
            lcoords = self.coords_to_crystal(lcoords, self.fmx, self.omx, self.n)
            zero_translation = {0:0, 1:13, 2:62, 3:171}[self.n]
        else:
            lcoords = array([lcoords])
            zero_translation = 0
        shape = lcoords.shape
        lcoords = lcoords.reshape((shape[0] * shape[1] * shape[2], shape[3]))
        kdt = kdtree(lcoords, bucket_size)

        # 0,0,0 translation is: Thickened cube numbers: a(n)=n*(n^2+(n-1)^2)+(n-1)*2*n*(n-1).
        # 1, 14, 63, 172, 365, 666, 1099, 1688, 2457, 3430, 4631, 6084, 7813 ...
        # python indexes -> -1
        asu_atoms = (shape[2] * shape[1])

        #QUERY
        result = defaultdict(dict)
        for idx, qcoord in enumerate(qcoords):
            ridxs, rdists =  kdt.rn(qcoord, search_limit)
            for (ridx, rdist) in izip(ridxs, rdists):
                entity_idx = ridx % shape[2]                    # 0 .. N -1
                symmetry = (ridx // shape[2]) % shape[1]        # 0 .. MXS -1
                translation = ridx // asu_atoms                 # 0 .. (2n + 1)^2 -1
                qent_id = qents_idx_to_id[idx]
                lent_id = lents_idx_to_id[entity_idx]
                if translation == zero_translation: # same asymmetric unit.
                    if mode == 'diff_asu':
                        continue
                    elif mode == 'diff_sym' and symmetry == 0:
                        continue
                    elif mode == 'diff_chain' and (qent_id[2] == lent_id[2]) and symmetry == 0:
                        continue
                result[qent_id][lent_id] = (sqrt(rdist), "%s_%s" % (translation, symmetry))
        return result